clear all; close all; clc;
t = 1023; M = 8*t + 3; 
m = 32;
N = 1e7;

y = nan(1, N); y(1) = 13;
for i = 2:N
y(i) = mod(M*y(i-1), 2^m);
end
x = y / 2^m;
figure(1)
hist(x);
xlabel('x');
xlim([-0.1 1.1]);
a = -1; b = 3; x = x*(b-a) + a;
dx = 0.01;
[h, xb] = hist(x, a-dx/2:dx:b+dx/2);
figure(2)
p = h / (sum(h) * dx);
plot(xb, p);
ylim([0 1.1]);
xlabel('x'); ylabel('p(x)');
